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1.  Introduction 

It  is  desired  to  find  the  level  of  spares  inventory  (y)  and  number 
of  repair  channels  (c)  to  support  an  operating  population  of  M items. 

An  operating  item  may  randomly  fail  and  require  repair.  Both  the  time  to 
failure  and  the  repair  time  of  any  item  are  independently  distributed  ex- 
ponential random  variables  with  parameters  X and  u , respectively. 

All  items  are  identical  and  independent. 

Generally  one  is  interested  in  knowing  how  many  spares  and  repair 
channels  (or  repairmen)  are  required  in  order  to  support  the  system  at  a 
certain  level  measured  in  terms  such  as  the  probability  that  the  popula- 
tion is  operating  at  full  strength  or  the  probability  that  a certain  per- 
centage of  the  population  is  operational  or  the  probability  that  the 
spares  pool  is  not  empty  upon  a request  for  a spare. 

The  classical  finite  source  queueing  theory  (also  known  as  the 
machine  repair  problems)  can  be  used  to  provide  the  steady  state  distri- 
bution of  the  number  of  units  in  and  awaiting  repair  which  is  required 
to  calculate  the  service  level  measures.  Further,  using  the  costs  of 
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spares  and  repair  channels,  a procedure  can  be  developed  to  find  the  op- 
timal combination  of  spares  and  repair  channels  to  meet  a specified  level 
of  service  (see  Reference  [6]). 

In  many  cases,  the  provisioning  for  spares  and  repair  capacities 
must  be  planned  in  advance.  Also,  in  many  new  systems,  operating  units 
are  added  over  an  initial  period  until  the  population  reaches  its  full 
strength.  An  example  of  this  is  a gas  turbine  engine  fleet  which  is  to 
build  to  a full  strength  of  256  ships,  starting  with  an  initial  fleet 
size  of  5 the  first  year  and  adding  ships  over  an  11-year  period.  When 
the  population  size  changes,  especially  if  new  units  have  different  mean 
times  to  failure  or  mean  repair  times  (which  is  likely  due  to  technologi- 
cal growth  and  learning),  transient  effects  might  be  significant.  In  Ref- 
erence 16],  it  was  assumed  that  the  population  reached  instantaneous  steady 
state  every  year  at  new  population  values.  We  concentrate  here  on  methods 
of  analyzing  transient  effects  so  that  situations  can  be  noted  for  which 
steady  state  assumptions  may  be  in  doubt  as  well  as  those  for  which 
steady  state  assumptions  may  be  adequate. 

2.  Calculation  of  Transient  Probabilities 

The  state  of  the  system  is  defined  by  the  number  of  items  in  and 
awaiting  repair,  which  can  vary  between  zero  and  M+y  . Letting  p„(t) 

represent  the  probability  of  the  system  being  in  state  j at  time  t given 
that  it  started  in  state  i (0  < i,j,  < my)  , tt  (t)  be  the  unconditional 

probability  of  the  system  being  in  state  i at  time  t,  P(t)  = (pj^t)} 
be  the  (Mfy+l)  x (my+1)  matrix  of  transition  probabilities  p^(t)  , 
ir(t)  be  the  M+y+1  component  vector  whose  elements  are  the  ir^(t)  , and 
tt(0)  be  the  initial  probability  vector  whose  components  tt^(0)  are  the 
probabilities  that  the  system  starts  in  state  i,  we  have 

TT(t)  = TT(0)P(t)  . (1) 

Denoting  the  steady  state  probability  vector  by  ir  , 

ir  = lim  ir(t)  . (2) 

t-*» 
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It  is  not  necessary  to  obtain  ir(t)  if  only  ir  is  of  interest  since 
this  can  be  determined  directly  by  steady  state  finite  source  queueing 
results  (for  example,  see  [5],  pp.  118-125).  However,  to  assess  transi- 
ent effects,  it  is  necessary  to  calculate  ir( t)  ; in  fact,  it  is  of  inter- 
est to  compare  ir(t)  to  ir  to  observe  the  speed  of  convergence  to  steady 
state. 

Solving  for  7r(t)  requires  solving  a set  of  Mty+1  differential 
equations  in  M+y  unknowns.  Obtaining  analytical  solutions  is  extremely 
difficult  and  in  most  cases  impossible  unless  M+y  is  very  small.  How- 
ever, a variety  of  procedures  are  available  which  can  yield  numerical  so- 
lutions to  the  differential  equations.  Reference  [4]  compares  several 
of  these  numerical  procedures  and  concludes  that  the  so-called  randomi- 
zation procedure  has  some  advantages.  This  procedure  is  presented  (al- 
though not  under  that  name)  in  Reference  [3]  on  page  46. 

For  the  purposes  of  study  in  this  paper,  the  randomization  pro- 
cedure appears  useful,  although  if  solutions  for  ir(t)  are  desired  for 
a large  number  of  t values,  the  Runge-Kutta  method  may  be  more  efficient. 
The  randomization  procedure  will  yield  P(t)  , from  which  ir(t)  can  be 
calculated  using  Equation  (1).  The  calculation  procedure  is  as  follows. 

Let  Q denote  the  intensity  matrix  (infinitesimal  generator)  of  the 
Markov  chain  of  this  birth-death  system,  that  is 


where 
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(M-i+y) A , (y  < i < y+M) 

0 , (i  > y+M) 


iy 

cm 


(0  < i < c) 
(i  > c) 


The  birth-death  differential  equations  then  can  be  written 

(t)  - TT(t)Q  , (3) 

where  it' (t)  is  the  vector  whose  components  are  ir^(t)  * dir^CO/dt  . 

The  randomization  procedure  gives  a method  of  calculating  P(t)  using 
Q . Denoting  the  i,jth  element  of  Q by  q^  , define  a scalar  a as 

a = max  q^  . 

Let  6 be  another  scalar  such  that 


6 > a . 


Then,  form  a matrix  P with  element 


PU 


as 


P - {PljJ  - I+)q, 

where  I is  an  (Mfy+1)  x (Mfy+1)  identity  matrix.  It  can  be  shown 
that  the  elements  of  P(t)  can  be  obtained  by 


Pij(t) 


-3t 


l 

n=0 


Bntn 
n!  rij 


(4) 


where  p^  is  the  (i,j)th  element  of  the  matrix  P(n)  and  P*n*  is  the 

^brix  P raised  to  the  power  n . Thus  the  computations  are  not  simple 
since  they  involve  raising  the  matrix  P to  the  nth  power  for  n-0,1,2,... 
Obviously,  the  numerical  procedure  must  truncate  n at  some  appropriate 
value  and  our  criterion  was  to  truncate  n at  a number  N , such  that 


where 


»«<«  - *'8‘  !04f  -if  • 


ir"  ' -in  n!  rij 
We  are  fortunate  here  in  being  able  to  guarantee  this  since 


pij(t) 


e-6c  j p<„> 

n=0 


n!  rij 


p>-8t 


V (Bt)  (n) 

u-HU  ">  « ' 


Now 


n=N+l  n!  ^ ’ n! 


n=N+l 


e-Bt  J p(»)  < .-Pt  £ . e-Bt  r Bt  _ £ (Bt^l 

L n=0  n>  J 

- 1 - e~Bt  l S&l  . 

rt  11  • 


n=0 


Thus 


^ \ n 

I 

n=0 

Hence  we  can  choose  the  smallest  N such  that  the  absolute  error 


p„<0  - p?.(t)  < i - *'BC  l ml  . 

iJ  „n  n! 


X - e-Bt  t ml  < c. 

n=0  n‘ 


(5 


Because  some  of  the  Ct)  may  be  very  small,  a somewhat  better  criter- 


ion might  be  to  choose  N such  that 


1 - e"Bt  f I§££ 
L-  n! 


n=0 


E-Bt  \ mi  w 


< e ; 


(6: 


n=0 


n!  rij 


that  is,  the  percentage  error  is  less  than  or  equal  to  e . This  is 
somewhat  more  complex  to  calculate  however. 
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This  procedure  was  programmed  on  an  IBM  370/148.  Analytical  solu- 
tions are  obtainable  for  the  special  case  M = y = c = 1 for  which  there 
are  only  three  states  (i=0,l,2).  The  three  differential  equations  in 
three  unknowns  given  by  Equation  (3)  can  be  solved  using  Laplace  trans- 
forms (see  Reference  [1]).  The  numerical  procedure  was  checked  against 
the  analytical  solution  for  the  special  case  and  found  to  be  quite  accu- 
rate. Table  I shows  the  comparison  using  criteria  from  Inequality  (5); 

_ 3 

that  is,  the  absolute  error  criterion  with  £ = 10  , and  10  . We  see 

TABLE  I 

COMPARISON  OF  NUMERICAL  AND  ANALYTICAL  SOLUTIONS 
M = y = c = 1;  A=  0.5,  y = 1.0 

Anal.  (Analytical);  Num.  (Numerical);  Numa:  e = 10  Num^:  e = 10  ^ 


rt 

X) 

o 

o 

rt 

v-/ 

P01(t> 

p02(t) 

Anal. 

0 1 

0 

0 

a 

Num 

1 

0 

0 

Num^ 

1 

0 

0 

Anal. 

x .7265 

.2227 

.0508 

Num3 

.7260 

.2224 

.0506 

Num'3 

.7266 

.2226 

.0508 

Anal. 

3 . 6010 

.2768 

.1222 

Num3 

.6003 

.2767 

.1222 

Anal. 

5 .5775 

.2839 

.1386 

a 

Num 

.5770 

.2837 

.1385 

Anal. 

io  -5715 

.2857 

.1428 

Num3 

.5710 

.2854 

.1427 

aRunning  time  to  generate  four  3*3  P(t) 
matrices  (t=l,3,5,10) : 15.78  sec. 

^Running  time  to  generate  one  3x3  P(t) 
matrix  (t=l):  3 min. 

- 6 - 
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that  the  numerical  procedure,  with  e = 10  , is  always  exact  to  the 

second  decimal  place  and  in  most  cases  to  the  third.  While  lowering  e 

to  10  ^ helps  somewhat  on  the  accuracy  (see  the  t=l  case  in  Table  I), 
the  added  expense  in  running  time  becomes  enormous  and  hence  it  was  felt 
-3 

that  10  was  satisfactory. 


3.  Results 

Figure  1 shows  the  results  of  a case  where  M = 10,  y = 4,  c=2, 

X = .15,  y = 1.0  . The  measure  of  system  service  used  is  the  probability 
that  the  entire  population  is  operational.  Denoting  this  by  A(t)  , we 
have 

y 

A(t)  = l tt  (t)  . 
i=0  1 

Two  cases  are  shown,  one  in  which  the  initial  conditions  had  all  machines 
"up,"  that  is,  no  machines  in  or  awaiting  repair,  one  machine  operating, 
and  one  in  the  spares  pool  [tt(0)  = (1,0,0)  ],  and  the  other  where  both  ma- 
chines were  "down,"  that  is,  one  in  repair  and  the  other  in  the  repair 
queue  [tt(0)  = (0,0,  1)  ].  We  can  see  that  in  the  first  case  steady  state  is 
approached  rather  quickly  while  in  the  second  the  approach  is  consider- 
ably slower. 


There  are  two  factors  influencing  the  speed  to  steady  state.  One 
is  the  initial  conditions,  that  is,  how  far  the  system  is  from  steady 
state  initially,  and  the  other  is  the  speed  at  which  the  transient  terms 
(terms  depending  on  t)  of  the  full  solution  to  the  differential  equations 
damp  out.  In  some  cases  it  is  known  what  the  initial  conditions  are  or 
are  likely  to  be,  but  in  many  cases  information  about  the  initial  state 
of  the  system  may  not  be  available  or  reliable.  In  order  to  divorce  the 
effect  of  the  initial  conditions  from  how  quickly  the  transient  terms 
damp  out,  we  decided  to  look  at  the  P(t)  matrix  itself.  We  became  in- 
terested in  two  facets  of  the  P(t)  matrix.  First,  we  wished  to  obtain 
information  about  the  behavior  of  the  p (t)  as  t increases  and  second, 

we  wished  to  determine  for  a given  t how  "far"  P(t)  was  from  its  steady 
state  value  P(°°)  , where 

- 7 - 
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P(oo) 


that  is,  a matrix  in  which  all  rows  are  the  steady-state  probability 
vector  it  . This  necessitated  the  development  and  study  of  single  value 
measures  for  the  "distance"  between  P(t)  and  P(°°)  . 


3.1  Behavior  of  p^Ct) 

We  know  P(0)  is  an  (M+y+1)  * (M+y+1)  identity  matrix  and  P(°°) 
is  the  (M+y+1)  x (M+y+1)  matrix  with  all  rows  equal  to  the  vector  tt  . 
We  are  interested  in  what  happens  in  between.  For  the  case  M = 10  , 
y = 4,c  = 2,A*  .05  , p = 1.0  , we  plotted  the  first  four  probabili- 
ties of  row  three  of  the  matrices  P(t)  , t*0,l,3,5,10  and  00  ; that  is, 
we  plotted  p2Q(t)  » P21(t)  * ^22^  * and  p23^  for  the  above  values 
of  t (see  Figure  2). 


It  can  be  seen  that  not  all  of  the  p^(t)'s  approach  their  steady 

state  value  monotonically.  The  diagonal  element  p22(t)  and  the  first 

element  p2Q  do,  p^Ct)  from  above  and  p2Q(t)  from  below.  But  the 

other  elements,  P2^(t)  and  P2g(t)  * start  at  zero,  overshoot  their 

steady  state  values,  and  then  appear  to  monotonically  decrease  toward 
them.  We  will  show  that  p^(t)  will  always  approach  from  above 


and,  although  it  has  not  been  proven,  we  conjecture  that  if 
TTj  , then  Pi;j(t)  > iTj  for  all  t > ^ ; that  is,  P^U)  , 
at  zero,  overshoots  it  at  most  once,  and  then  approaches 

above . 


Pij 


(tx)  = 

, starts 
from 


The  above  conjecture  can  be  argued  by  the  following  physical  anal- 
ogy. Consider  a row  of  the  P(t)  matrix  as  a bar  of  water.  At  t-0  , a 
sugar  cube  is  dropped  in  at  the  diagonal  (i,i)  position.  As  time  passes. 


- 9 - 
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the  sugar  is  dispersed  along  the  row  until  it  reaches  its  steady-state 
concentration  (i=0,l, 2, . . . ,FHy)  . Note  that  at  i,i  the  sugar 

gradually  flows  away  from  its  initial  concentration  (of  1)  to  its  final 
concentration  . At  either  end,  sugar  gradually  flows  in  until  it  builds 

up  to  the  steady  state  concentrations  irrt  and  it...  . All  other  states 

U M+y 

start  at  zero,  receive  sugar  from  the  diagonal  element,  some  of  which  will 
eventually  flow  past  to  other  states  and  some  of  which  will  remain,  causing 
the  overshoot  prior  to  settling  down  to  its  steady-state  concentration. 

It  is  possible  that  other  states  near  to  the  "end"  of  the  row  may  not 
overshoot  depending  on  the  rates  of  sugar  flow  in  versus  sugar  flow  out. 

We  expect  the  greatest  overshoots  to  be  in  the  states  near  the  diagonal 
element. 

We  next  looked  at  the  behavior  of  the  p^Ct)  38  a function  of  j 

for  fixed  t . That  is,  it  was  of  interest  to  observ-  which  p„(t)  were 

larger  than  their  steady  state  values,  ir^  , and  which  were  smaller,  for  a 

given  value  of  t . In  all  cases  run,  we  observed  one  of  three  situations 
as  shown  in  Figure  3.  It  can  be  seen  that  the  p^CtJ's  which  overestimate 

their  tt^'s  appear  together  in  a "clump."  If  this  were  always  true,  there 

could  be  significant  savings  in  the  computational  effort  required  for  ob- 
taining the  distance  measures  ( | P(«)— P(t ) | ) discussed  in  the  next  section. 

We  have  proved  that  this  clumping  will  always  be  the  case  as  well  as  its 
relative  position  in  the  row  in  the  lemma  and  theorems  presented  below. 

Lemma:  Let  N^(t)  * the  number  of  customers  in  and  awaiting  re- 
pair at  time  t given  i in  and  awaiting  repair  at  time  0.  Then 

Pr(Ni+1(t)  < k)  < Prj^Ct)  < k)  < Pr(N1+1<t)  < k+l)  . 

Proof:  Let 

xi  = min{u:  N^u)  * Nj_+i^u)  » 0 < u _<  00 } . 

It  is  obvious  that 
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Pr^CO-k  | J±<t]  = Pr[Ni+1(t)=k  | T±<t]  , k=0,l,2,...  , (7) 

Pr[Ni(t)<k  | T±>t]  > Pr[N1+1(t)<k  | T±>t]  . (8) 

Combining  (7),  (8),  and  the  fact  that  Pr(xi>t)  > 0 yields 

Pr(Vt)ik)  > Pr(Ni+1(t)  < k)  . (9) 

Now  we  denote  by  A^(t)  the  number  of  arrivals  in  [0,t]  and  by  S^t) 

the  number  of  service  completions  in  the  same  time  interval,  given  i in 
system  at  t=0  . Since  > AJ+1  and  < y for  all  j > 0 , then 

PrjA^Ct)  - S^t)  < k)  < Pr(Ai+^(t)  - si+-^(t)  ,<  k ) , for  any  k . 

The  proof  is  now  complete  since 

Nj (t)  = j + A^t)  - (t)  . 

From  (9)  we  obtain  that 

P00(t)  = pio^fc^  = p20^  = *'*  * 

and  as 

M+Y 

*0  ■ ±l  Vio<‘>  • 

then 

p00(t>  i ”0  • 

The  next  theorem  extends  this  result  for  all  Pli(t)  . 


Theorem  1 : 


pii(t)  > \ 


v i and  t. 


Proof:  We  compare  our  system  to  a Markovian  queueing  system 

with  only  two  states,  0 and  1.  Denote  the  rate  out  of  state  0 by  V and 
the  rate  out  of  state  1 by  n • Let  rQ  ±(t)  be  the  transient  probability 

of  being  in  state  i (i=0,l)  at  t given  that  the  system  started  in  state  0. 
The  steady  state  probabilities  of  the  two-state  system  depend  on  the  ratio 
v/n  . We  choose  v/n  such  that 
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lira  rQ  (t)  = ir  . 

t-*»  U 1 

It  is  well  known  that  rQQ(t)  > t\±  for  any  V , n , and  t ; hence  the 
required  result  is  obtained  if  n and  v can  be  chosen  to  yield 

rQ0(t)  < Pii(t)  . (1C 

Since  rQ()(0)  = p^O)  = 1 , then  Equation  (10)  is  satisfied  if 

r00(t)  = pii(t)  d* 


at  any  point  t where 


r00^  pii^t^  * 


To  determine  the  conditions  under  which  Inequality  (11)  holds,  we  write 
Pii(t)  = PIi(t)  - pIi(oo)  = -(^i'Hi1)(Pli(t)-iri) 

(13) 

and 

r00(t)  " r00(t)  • r00(oo)  = -(^)(roo(t)-7ri)  • (14) 

Using  Equation  (12)  in  (13)  and  then  comparing  (13)  to  (14)  yields  that 
Inequality  (11)  holds  iff 


-(\hK))  < 


-Ul*Ul>|P11(t)-i'1|  + >1.1lP1.1.1(t)-ir<.1|  + U1+1|P,  I+1(«-Vl) 


pll(t)  - "i 


The  RHS  of  (15)  is  clearly  bounded  since  the  matrix  P(t)  can  be  ex- 
pressed as 

V Cit  C^t 

P(t)  - lira  P(t)  + B e 1 + B,e  3 + ...  + B_1.ve^!+Y  , 

t-xx>  z J 

where  < 0 is  the  ith  eigenvalue  of  the  infinitesimal  generator 
matrix  Q and  Bi  is  a (M+Y+l)  x (MfY+1)  matrix  that  has  finite 
elements  (see,  for  example,  [2],  pp.  364-368).  Thus  we  conclude  that  V 
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and  n can  be  chosen  big  enough  so  that  (11)  is  satisfied  at  any  point 
at  which  (12)  is  satisfied,  and  the  proof  is  now  complete. 


Theorem  2:  For  any  i and  t there  exists  a j^(i,t)  and  a jf(i,t)  , 


j^i.t)  < i < Jr(i,t)  , such  that 


Pij(t)  * TTj  > 0 , if  j * j£(i»t)  or  j = jr(i,t) 

Pi:j(t)  - ir  > 0 , if  j£(i,t)  < j < jr(i,t) 


p^(t)  - tt^  < 0 , otherwise. 


Proof:  Denote 

K(i,t)  = (k:  P±k(t)  = TTfc}  . 

At  t=0  we  have  p^^O)  “ 1 and  p^(0)  = 0 , i^j  ; hence  the  required 

conditions  are  satisfied.  To  complete  the  proof  we  have  to  show  that  for 
any  k e K(i,t)  we  have 

Pik(t)  - 0 * if  k < jA(i.O-l  or  k > jr(i,t)+l 

Pik(t)  - 0 ’ if  j£(1>t)  < k < jr(i,t)  * 

To  show  that  the  above  conditions  hold,  we  write 

pik<c)  * pik(t)  * pIk(">  ' 


M+Y-l  > 


-Xk^Pik(t)"\J  + ^k+l l Pi , k+1  ^ -7Tk+l ^ * k=° 

'(Ak+\^Fik(t)"\J  + Xk-l^Pi,k-l(t)~\-l^  + \+llpi,k+l(t)-Trk+l)  ’ k > 0 

-\<pik(t)-\)  + VilPi.k-^^-Vi*  > k=MfY 


, k=M+Y 


If  k < j^(i,t)-l  or  k > jr(i,t)+l  , then  p±  kl(t)  - ir^  < 0 and 

pi,k+l(t^  " \+i  = 0 » thus  pik^^  = 0 * If  < k < JyCitt)  » 

then  P^.^O  - Vl  > 0 and  Pi>k+1(t)  - *k+l  > 0 ; thus  P^(t)  > 0 . 
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3.2  Measures  of  the  distance  between 
P(t)  and  P(°°) 


We  now  turn  our  interest  to  gaining  information  concerning  how 
far  from  steady  state  the  system  is  at  a given  time.  We  desire  some 
measure  of  the  distance  between  two  matrices  P(t)  and  P(»)  ; that  is, 
a measure  of  |P(t)  - P(°°)  | . 


Consider  the  following  example  for 


P(tx)  , P(t2) 


and  P(°°)  : 


velop  scalar  measures  to  reflect  this  situation. 


Two  measures  that  appear  reasonable  are  the  following: 

Mx(t)  = M1(t)/M1(0)  , 

where 

M+y  M+y 

V0  = l TT  l |p  (t)  - TT  | (17) 

i-0  1 j=0  J 

and 

M2(t)  = M2(t)/M2(0)  , 

where 

M+y  M+y 

Mo(t)  "ll  • (18) 

i=0  j=0  1J  j 


M^(°°)  - 0].  The  second  is  merely  the  sum  of  all  the  absolute  differ- 
ences between  the  corresponding  elements  of  P(t)  and  P(°°)  . The 
first  is  a weighted  sum  of  the  corresponding  element  absolute  differ- 
ences, the  weights  being  the  steady  state  probabilities;  that  is. 


i 
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errors  in  the  state  probabilities  which  are  seldom  visited  in  steady 
state  are  "discounted." 


A well-known  mathematical  definition  of  the  distance  between  two 
matrices  is  the  norm  of  the  matrix  formed  by  subtracting  corresponding 
elements,  denoted  by  | |P(t)  - P(<»)  | | . Matrix  norms  must  obey  certain 
properties  (see  Reference  [7],  page  59),  which  are  given  below. 

(1)  | |P(t)  - P(»)  | | > 0 and  | |P(t)  - P(<*>)  ||=0, 


iff  P(t) 

- P(oo) 

- o ; 

(2) 

| |a[P(t) 

- P(»)J 

II  - |a|  • 

l|P(t) 

- P(»)  | | ; 

(3) 

lltP(t1) 

- P(“)l 

+ [P(t2) 

- P(»)  1 | 

| < | |P(t1)  - P(»)  | | + 

l|P(t2)  - 

P(«)  1 1 

9 

(4) 

- P(“)l 

• tP(t2)  - 

P(«)  ] | | 

< | |P(tx)  - P(«)  | | • 

1 lp(t2)  - 

P(co)  | | 

• 

We  have  not  been  able  to  show  that  and  M2  obey  all  the  above 

properties,  although  we  believe  that  they  do.  Howevei , variations  of 
and  M2  which  do  obey  all  properties  are  (see  Reference  [7], 

page  61): 


M3(t)  = M3(t)/M3(0) 

where 


M+y 

M,(t)  = l tt  max  |p . . (t)  - tt  | (19) 

i=0  1 j 1J  3 

and 

M4(t)  = M4(t)/M4(0) 

where 


M+y 

M (t)  = l max  |p.,(t)  - it  | . (20) 

j=0  i 13  3 

Once  again,  M3  uses  the  rr^  to  weight  the  maximum  row  deviations, 
while  M4  is  based  on  the  unweighted  maximum  column  deviations.  We 
have  run  many  cases  to  observe  the  behavior  of  the  proposed  measures. 
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We  chose  situations  for  which  M+y  equaled  five  so  that  we  could  view 
the  entire  6x8  P(t)  matrix  and  compare  it  with  the  6x6  P(<»)  matrix. 
The  cases  run  are  shown  in  Table  II,  with  the  results  shown  in  Table  III. 


TABLE  II 
6x6  CASES  RUN 


Before  comparing  the  four  measures  numerically,  we  would  like 
to  present  two  properties  that  may  help  in  judging  the  measures. 


WfY 

Property  1:  £ |p  (t)  - it  | 

j=0  3 3 

This  follows  from  Theorem  2,  which  yields 


2 l (Pi1(t)  - n ) . 
j=j.(i,t)  iJ  3 


M+Y  jr(i,t) 

l |p44(t)-irj  * l (p±JCt)  - ttJ  + 


j-o 


ij 


j 


j 


l K - Pij(t)l 


j<j£(i,t) 


+ l 


( TTj  - P^tt))  ; 
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TABLE  III 


RESULTS 

OF  6x6 

CASES 

RUN 

A(°°) 

Measures 

t-0 

t-1 

t*»3 

t-5 

t-10 

.820 

*1 

1.000 

.376 

.143 

.066 

.010 

M2 

1.000 

.505 

.232 

.115 

.019 

M3 

1.000 

.319 

.103 

.046 

.007 

*4 

1.000 

.366 

.178 

.092 

.016 

.820 

“I 

1.000 

.216 

.046 

.010 

.001 

M2 

1.000 

.335 

.081 

.019 

.001 

M3 

1.000 

.170 

.031 

.007 

.000 

1.000 

.246 

.066 

.016 

.000 

.820 

*1 

1.000 

.142 

.015 

.002 

.001 

M2 

1.000 

.232 

.027 

.003 

.000 

M3 

1.000 

.103 

.010 

.001 

.000 

«4 

1.000 

.178 

.023 

.003 

.000 

.820 

“l 

1.000 

.097 

.005 

.001 

.001 

M2 

1.000 

.164 

.009 

.001 

.001 

M3 

1.000 

.068 

.003 

.000 

.000 

5« 

1.000 

.129 

.008 

.000 

.000 

.964 

*1 

1.000 

.266 

.039 

.007 

.001 

M2 

1.000 

.464 

.110 

.025 

.001 

M3 

1.000 

.262 

.039 

.007 

.001 

*4 

1.000 

.325 

.092 

.022 

.000 

.630 

\ 

1.000 

.221 

.043 

.009 

.001 

M2 

1.000 

.275 

.056 

.012 

.001 

*3 

.020 

m 

Ma 

.099 

■ ■-V  ■■  - : ■ 
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Case 

A(°«>) 

Measures 

t-0 

ww 

Ei 

t=10 

7 

.462 

*1 

1.000 

.198 

.033 

.006 

.000 

*2 

1.000 

.237 

.041 

.007 

.000 

53 

1.000 

.103 

.016 

.003 

.000 

5, 

1.000 

.164 

.023 

.004 

.000 

8 

.624 

«1 

1.000 

.445 

.223 

.123 

.034 

M2 

1.000 

.515 

.276 

.159 

.044 

«3 

1.000 

.371 

.120 

.062 

.016 

1.000 

.417 

.199 

.111 

.031 

9 

.462 

“l 

1.000 

.478 

.258 

.154 

.052 

«2 

1.000 

.536 

. 30 '♦ 

.186 

.063 

»3 

1.000 

.416 

.145 

.077 

.025 

1.000 

.457 

.215 

.125 

.037 

10 

.860 

\ 

1.000 

.371 

.155 

.076 

.013 

M2 

1.000 

.522 

.268 

.141 

.026 

»3 

1.000 

.332 

.128 

.061 

.011 

1.000 

.384 

.205 

.114 

.022 

11 

.735 

“l 

1.000 

.349 

.121 

.053 

.007 

M2 

1.000 

.535 

.260 

.127 

.019 

*3 

1.000 

.319 

.105 

.044 

.006 

«4 

1.000 

.400 

.211 

.111 

.017 

12 

.656 

*L 

1.000 

.359 

.114 

.048 

.006 

52 

1.000 

.521 

.227 

.106 

.014 

*3 

1.000 

.315 

.092 

.036 

.004 

5< 

1.000 

.376 

.186 

.092 

.013 
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and  as 


1 - l P^Ct)  - I ir  , 
j J j J 


rrri 

L |p±1(t)  - ir  | 
j=0  ^ J 


2 l (P±1  (t)  ‘ "J 


Property  1 suggests  that  M^(t)  and  M2(t)  may  in  general  re- 
quire less  computation,  since  only  part  of  the  matrix  P(t)  has  to  be 
calculated.  Further,  while  all  numerical  work  indicated  that  j£(i,t) 

and  j (l,t)  do  not  decrease  as  i increases  (that  is,  the  "clump," 
P^jC^)  ” 77 j = > moves  to  the  right  as  we  go  down  rows  of  the  P(t) 
matrix),  we  have  not  been  able  to  prove  this. 

Property  2:  The  measures  M^t)  , M2(t)  , and  M4<t)  decrease 
monotonically  in  t. 

To  show  that  M^(t)  and  M2(t)  are  monotone,  we  use  Expression 
(16)  to  obtain  for  j£(i,t)  t 0 and  Jr(i,t)  t M4Y  that 

* j-jjl.t)  |Pij(t)  ■ ^ = ,(i,t)-i(t)  - ^j^i.o-i) 

\<i,t)(Pi,j£(i,t)(t)  ' \(i,t)) 

Ajr(i.t)(Pl,jr(i,t)(t)  " 7rjr(i,t)J 

+ Ujr(i,t)+l(Pi,jr(i,t)fl(t)  " Trjr(i,t)+l) 

(21) 

Expression  (21)  yields  the  monotone  convergence  since  all  four  terms  of 
the  RHS  are  nonpositive.  The  results  for  the  cases  j£(i,t)  = 0 and 

Jr(i»t)  * M+Y  can  be  obtained  in  a similar  way. 
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We  now  show  that  M^,(t)  is  monotone.  First  we  notice  that 
(t)  | = max  maxtp^t)}  - n.  ; -jminlp^Ct) } - Tr^jj  . 

Thus  the  statement  is  true  if  max(p  . ( t) } is  decreasing  and 

i 1J 

(t)}  is  increasing.  These  requirements  are  met  since 
i 

p(tl+t2)  = P(t1)P(t2)  ; 

hence  the  elements  of  the  jth  column  of  P(t^+t2)  are  convex  combina- 
tions of  the  jth  column  of  P(t2>  . 

As  for  M^(t)  , in  all  our  calculations  this  measure  was  monotone 
decreasing.  However,  we  are  unable  to  prove  it. 

Since  we  are  interested  in  assessing  the  relative  merits  of  all 
the  measures,  fjr  the  cases  in  Table  II  all  p^(t)  were  calculated. 

Case  1 can  be  considered  the  base  case.  Cases  2,  3,  4 increase 
A and  y , keeping  p (defined  as  MA/cy)  constant.  Cases  5,  6,  and 
7 have  y constant  but  varying  A , while  cases  8,  9,  and  10  have  A 
constant  with  varying  y . Cases  11  and  12  revert  back  to  the  base 
case  p , but  vary  M , y , and  c . 

Figure  4 shows  some  graphs  of  the  measures  as  they  vary  with  cer- 
tain parameters.  In  4(a),  we  see  that  all  measures  do  decrease  monoton- 
ically  with  time  (in  fact,  they  appear  convex),  which  certainly  must  be 
a requirement  for  the  measures  to  be  useful.  The  unweighted  measure  M2 

decreases  more  slowly  than  its  weighted  counterpart,  . The  same  is 
true  of  and  . Of  course,  is  always  less  than  , and 

is  always  less  than  M2  , since  and  consider  only  maximum  row 

and  column  errors,  respectively. 
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In  Figure  4(b),  we  see  Chat  increasing  X anrl  U together  so  as 
to  keep  p constant  causes  a faster  approach  to  steady  state  (all  mea- 
sures again  decrease  monotonically) . This  is  intuitive , since  we  would 
expect  the  greater  the  transition  (arrival  and  service  completion)  rate, 
the  faster  steady  state  should  be  approached.  While  this  generally 
holds  also  for  varying  only  arrival  rate  (service  completion  rate) 
keeping  service  completion  rate  (arrival  rate)  fixed  [see  Figures  4(c) 
and  4(d)],  the  measures  do  not  always  decrease  monotonically.  In  the 
case  where  y is  fixed,  decreases,  then  slightly  increases,  then 

decreases  again  as  X increases.  In  other  runs  not  shown  in  Table  III  we 
saw  similar  effects  for  the  other  measures.  Just  why  this  is  so  is  not 
intuitively  clear,  although  one  must  keep  in  mind  that  the  traffic  in- 
tensity (p)  is  now  also  varying,  so  that  the  probability  distributions 
of  the  system  states  are  also  changing. 


The  final  two  cases  run  (Cases  11  and  12)  had  the  same  p but 
varied  M , y , and  c . Even  with  p fixed,  A(«)  varies,  showing 
the  effect  of  population  size,  spares  pool  size,  and  repair  capacity 
on  availability.  However,  the  convergence  to  steady  state  as  given  by 
the  four  measures  does  not  appear  to  change  significantly. 

Which  measure  to  choose  and  how  to  utilize  the  measure  ultimate- 
ly lie  in  a subjective  judgment.  We  see  that  M2  decreases  more  slowly 

than  the  others,  while  decreases  most  rapidly.  Since  in  effect  the 

scale  is  arbitrary,  one  a measure  (or  measures)  is  chosen,  subjective 
benchmarks  would  have  to  be  set  up.  For  example,  consider  the  P(t) 
matrices  for  Case  1,  t=l,3,5  and  <*>  , and  the  respective  measures  shown 
in  Figure  5.  The  full  P(t)  matrices  for  t=*l,3,5  are  given.  At  the 
bottom  of  each  is  listed  the  steady  state  probability,  [the  P(°°) 

matrix  has  all  its  rows  identical  and  equal  to  it  « (it  ,ir  , . . . ,tt„.v) 

0 1 rrrY 

so  that  near  steady  state  all  elements  in  column  j of  a P(t)  matrix 
should  be  close  to  tt^]. 
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Figure  S. — P(t)  , M (t)  and  A(t)  for  Case 
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In  viewing  Figure  5,  it  is  fairly  clear  that  P(l)  is  extremely 
far  from  P(°°)  and  P(3)  is  quite  far,  while  P(5)  is  somewhat  closer — 
although  in  a few  places  there  are  some  sizable  errors,  particularly  in 
rows  5 and  6,  which,  of  course,  are  discounted  by  measures  and  . 

Thus  one  might  conclude  that  a value  of  of  .05  or  less  indicates  a 

near  steady  state  condition,  while  a value  for  of  .10  or  less 

suffices. 

Also  shown  in  Figure  5 is  the  steady  state  population  availability, 

A(°°)  , along  with  two  transient  A(t)'s  . One  of  these,  ACt)^  » as~ 

sumes  all  units  are  operational  at  time  zero,  and  the  other,  ACt)^  , 

assumes  the  opposite,  namely,  all  units  are  down  (in  or  awaiting  repair) . 

Note,  as  before,  that  A(t)^p  is  much  closer  to  A(°°)  than  is  A(t)DN  , 

since  A(»)  is  nearer  to  A(0)im  *=  1 than  to  A(0)„„  = 0 . 

Ur  DN 

As  to  which  measure  might  be  preferable,  again  subjectivity  is 
required.  Consider  Figure  6,  where  P(l)  for  Cases  1,  11,  and  12  is 
shown,  along  with  the  respective  P(°°)  , A(t)  , and  A(°°)  . In  viewing 
the  matrices,  it  appears  that  all  are  about  the  same  distance  from 
steady  state.  In  looking  at  A(t)^,  and  A(t>DN  versus  A(<*>)  , perhaps 

one  might  say  that  Case  12  is  the  "poorest”  case,  while  Case  1 is  the 
"best."  None  of  the  measures  indicates  Case  12  to  be  the  poorest,  al- 
though M2  and  come  closest  to  this.  The  measures  and 

do  indicate  that  Case  1 is  the  best. 

1 

In  Figure  7 two  more  cases  are  compared,  namely.  Cases  2 and  6. 

Again  it  is  difficult  to  judge  which  P(l)  is  the  better,  although  in 
looking  at  the  A(l)'s  one  might  give  "the  edge"  to  Case  2.  Here, 
however,  shows  Case  2 to  be  preferable,  while  M2  , , and 

show  the  opposite. 

It  appears  that  no  measure  is  superior,  so  that  the  choice  should 
perhaps  be  the  one  that  is  computationally  easiest  to  obtain.  In  view  of 
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the  theorems,  and  M2  are  computationally  most  efficient  for  large 

matrices,  so  that  either  of  these  might  be  an  appropriate  choice. 
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